# In this script we price a bond according to the market yield and estimate its
# Present Value Sensitivity using Duration and Convexity using continuous 
# compounding. 



library(tidyverse)

# Market yield (computed in`S3.Estimating the market yield curve by replication.R`)
load("Examples/market.yield.curve.Rdata")
market.yield.curve

# Let's use continuous compounding.
market.yield.curve <- market.yield.curve %>% select(t, P_t, y_t = y_t.cont)

# Let's consider the following bond.
coupon.rate <- 0.05
face.value <- 100
maturity <- 5

# Price the bond (note that I add a first line full of zeros)
bond.market.price.tbl <- market.yield.curve %>% 
  head(n = maturity) %>% 
  bind_rows(data.frame(t = 0, P_t = 0, y_t = 0)) %>%
  arrange(t) %>% 
  mutate(
    # Payments of our bond
    c_t = c(0, rep(coupon.rate*face.value, maturity)) + c(rep(0,maturity), face.value),
    # The bond price using the market yield
    `PV of c_t` = P_t*c_t
  )
bond.market.price.tbl

save(bond.market.price.tbl, file = "Examples/bond.market.price.tbl.Rdata")

# Computing the sensitivity measures - Duration and Convexity
sensitivity <- bond.market.price.tbl %>% 
  mutate(
    # In accordance to the formulas in the slides 
    d_aux = t * P_t * c_t,  # or else, t * exp(-y_t * t) * c_t
    c_aux = t^2 * P_t * c_t # or else, t^2 * exp(-y_t * t) * c_t
  ) %>% 
  summarize(
    bond.price = sum(`PV of c_t`), d_aux = sum(d_aux), c_aux = sum(c_aux)
    ) %>% 
  mutate(duration = d_aux/bond.price, convexity = c_aux/bond.price)
sensitivity

# Now let us consider a parallel shift of 1% (\delta\bar{y} = 0.01) and 
# approximate the bond price accordingly.
paral.shift <- -0.01

sensitivity <- bond.market.price.tbl %>% 
  mutate(
    d_aux = t * exp(-y_t * t) * c_t,
    c_aux = t^2 * exp(-y_t * t) * c_t
  ) %>% 
  summarize(bond.price = sum(c_t*P_t), d_aux = sum(d_aux), c_aux = sum(c_aux)) %>% 
  mutate(duration = d_aux / bond.price, convexity = c_aux / bond.price)


sensitivity

# Now let us consider a parallel shift of 1% (\delta\bar{y} = 0.01) and 
# approximate the bond price accordingly to first and second order.
paral.shift <- -0.01

sensitivity <- sensitivity %>% 
  mutate(
    approx.1st.order = bond.price - bond.price * duration * paral.shift,
    approx.2nd.order = approx.1st.order + .5 * convexity * paral.shift^2
    )

sensitivity$approx.1st.order
sensitivity$approx.2nd.order

# Let us check what would be the exact new value (homework - similar to what we
# have done in script '2.Example bond valuation.R')










pert <- c(0,rep(1,5))*paral.shift
pert
bond.price.perturbed.1 <- bond.market.price.tbl %>% 
  mutate(
    y_t = y_t + pert,
    P_t = exp(-t*y_t), # remember we are using continuous compounding
    `PV of c_t` = P_t*c_t
    )
bond.price.perturbed.1
bond.price.perturbed.1 %>% summarise(sum(`PV of c_t`))

# Repeat the same for a parallel shift of 1% (homework too).
paral.shift <- 0.01

approx.1st.order <- sensitivity$bond.price * (1 - sensitivity$duration * paral.shift) 
approx.1st.order

approx.2nd.order <- sensitivity$bond.price * 
  (1 - sensitivity$duration * paral.shift + .5 * sensitivity$convexity * paral.shift^2)
approx.2nd.order


pert <- c(0,rep(1,5))*paral.shift
pert
bond.price.perturbed.1 <- bond.market.price.tbl %>% 
  mutate(
    y_t = y_t + pert,
    P_t = exp(-t*y_t),
    `PV of c_t` = P_t*c_t
  )
bond.price.perturbed.1
bond.price.perturbed.1 %>% summarise(sum(`PV of c_t`))